The Extended Plane Wave Expansion Method in Three Dimensional Anisotropic 

Photonic Crystal 
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In this paper, we extend the conventional plane wave expansion method in 3D anisotropic photonic 
crystal to be able to calculate the complex k even if permittivity and permeability are complex 
numbers or the functions of u. There are some tricks in the derivation process, so we show the process 
in detail. Besides, we also provide an example for testing and explaining, and we also compare the 
results with the band structure derived from conventional plane wave expansion method, then we 
finally find that there is a good consistency between them. 
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Recently, the researches of the properties of the 
photonic crystals (PCs) have aroused great interests, 
since the concept of the PCs has been proposed by 
Yablonovitch and JohnQ, 0, 0- Briefly speaking, PCs 
are periodically structured electromagnetic media, gen- 
erally processing photonic band gap (PBG). Most of the 
studies stress the PBG structures with the use of conven- 
tional plane- wave expanded (PWE) method0,0|. How- 
ever, there are still many articles explore the influence of 
interface, such as the studies of transmission, reflection, 
and the penetration depth etc. |E IS Furthermore, 
the penetration depth relates to the imaginary part of 
wave vector. As for the complex k calculation in 2D 
isotropic photonic crystals, we had sufficiently discussed 
about it in the last paper[13]. Now, this paper is to con- 
tinue with the last one. Furthermore, the emphasis of 
this paper is put on the general formula, 3D anisotropic 
case, of extended plane wave expansion (EPWE) method. 

Though the main part of the idea resembles in 2D 
isotropic case [13], the formula and derivative process are 
much more complicated than that in 2D isotropic case, 
because the basis of wave functions can not be treated as 
scalar functions, TE and TM modes in 2D isotropic case. 
However, the problem of the difficult part has been over- 
come and we will explain it in the following description. 
Besides, the eigenfunctions set derived from this EPWE 
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method is completely the same as that derived from the 
conventional PWE method. So we have no qualms about 
the inaccuracy of the propagation modes between these 
two methods. 

The system we discussed is periodically structured 
without charge p and current J. Therefore, according to 
Maxwell Equation, the magnetic field H(r) should obey 



(k + G) x e^_ G , (k + G') x He = w 2 /i G _ G <H G , 



(1) 
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G and G' are the reciprocal lattice vectors, u> and k 
are the frequency and wave vector, e(r) and fi(r) are the 
tensors of permittivity and permeability of which e G and 
/t G are the Fourier expansion components, respectively. 

Now, let us expand Eq.(l) directly through x,y and z 
directions 



x: 
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y: 



(k + G) z e-l (k + G') z + (k + G) x e" 1 (k + G') y - (k + G) z e" 1 (k + G% - (k + G) x (k + G') z \ H x 
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(2b) 



(* + G )x tyy ( fc + G ')z + ( fc + G ) y (* + G\ - (k + G) x e" 1 (k + G') y - (k + G) y (k + G') 
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= — W (fl zx H x + ^zyHy + fi zz H z ) , 
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where and fiij are the abbreviations of £g-g'. i,j and 
AG-G',i,j ; an d -ffi is the abbreviation of i?G',i- When 
k is provided, Eq.(2) becomes an eigenvalue problem in 
which the eigenvalue is uj and is the conventional PWE 
method. Now, there comes up an interesting question 
that is whether k must be a vector of which the com- 
ponents are real numbers. The answer is "No", and we 
just need to do some modification on Eq.(2) to get the 
complex k, because Eq.(2) is a 4 variables (k and u) 
equation. 

In the beginning, two important things need dis- 
cussing. First, the inner product of (k + G) and Eq.(2) 
results in 



(h + Gi)^G-G',i,jH G ,j 







G',i,j 



which are the restriction functions of which the amount 
is N, meanwhile, N is the amount of {G} set. Therefore, 
the certain amount of the independent eigenfunctions in 
Eq.(2) is 2N not 3N. That's why we will get the fake 
eigenvalues which are uo 2 — if Eq.(2) is calculated as an 
eigenvalue equation directly. 

To avoid this situation occurring in our method, the 

eigenvector we selected in our method is ' 



H 



GX 



not 



H G 
H G 



where Hqi and Hq arc 



H G,y 
Hg.z 



and 



Hg,x 



, Hgi and Hq arc ^Hgi and fc^Hc, re- 



Hgi 

spectively. 

Second, there are no k x HG',% , i = x,y, z, and k x HG',x 
in Eq.(2a), which is the x component of Eq.(l), because 
the inner products of Eq.(l) and x will cause the exis- 
tence of just one k x or even no, and (k + G/) x Ho part 
will restrict the existence of k x HG, x - 



Therefore, the treatment of x component will be dif- 
ferent from y and z components. The following is the 
detail derivation process: 

First of all, the y and z components of Eq.(2) can be 
written as a matrix formula 



" - . ~ . ~ . . 1 / H G 
Bi:B 2 :Ci:C 2 I jj 

and its expansion type is 
Bii?G,rr + CiHg,x + B 2 :C 2 



= A (fc x H G x) 



Hgi 
Hgi 



k x k (h G x) 

(3) 



where A, Bi, B 2 , Ci, C 2 are 2Nx2N, 2NxN, 2Nx2N, 
2N x N and 2N x 2N matrices and their elements will 
be illustrated in Appendix. 

As regards the x component of Eq.(2), we can write in 
another form which is different from y and z components 
of Eq.(2). Thus the matrix form of Eq.(2a) is 



EiiEa) (H g ) = d(h g ±) 



(4) 



where D, Ei, E 2 are Nx2N,NxN,Nx2N and their 
elements are also in Appendix. 
From Eq.(4) we obtain 

H G ,x = -Er 1 E 2 H G ± + VDHgx, (5a) 
H G , X = -Ej-^aHcx + fc.E^DHcx, (5b) 

where Eq.(5b) is the production of Eq.(5a) multiplied by 
k x . A combination of Eqs.(3) and (5) yields 
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A-CiEi'D 
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BiEr 1 
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where F is a 2N x AN matrix. Considering the equation 
above with fcjHci = Hgi, wc finally have an equation 




which is an k x eigenvalue equation, and the order of 
eigenfunction ( ) is 4N. In addition, and I are 

2N x 2N zero matrix and identity matrix, alternatively. 
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FIG. 1: The schematic view of a cubic unit cell in which three 
GaAs square rods cross together from the x, y and z direction. 
The lattice constant, width of square rods and e of GaAs are 
a, 0.4a and 11.43eo, respectively. 



find that lu = 0.2-^p is not located in band gap, so such 
kind of condition should also appeared in our method 
when we choose the same lu to plot the contour line or 
surface. Figure (2b) in which lu = 0.2^l£ , k z = and k y 
scanned from — — to - is the figure of real value solution 

a a ° 

of k x derived from Eq.(6). When Fig. (2a) compares with 
Fig. (2b), we will find out the width of contour in Fig. (2b) 
equals the width of X — ■> T region in Fig. (2a). 
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For testing this method, we use an Intel centrino 1.4G, 
512 MB RAM with matlab code published on mathworks 
website to run an isotropic simple cubic case in which 
the GaAs square rods — their widths are 0.4a, and a 
is the lattice constant — cross together from x,y and z 
direction in the vacuum. In this system the permittivity 
e of GaAs and vacuum are 11.43eo and eo, alternatively, 
and the permeability is /io everywhere. You can see 
its structure in Fig.(l) and calculation results in Fig. (2). 
We spent about 6 hours on getting Figs. (2b) and (2c) 
when using 729 {G} and taking 17 k y points from to 
— to accomplish the calculation. As regards Fig. (2a), it 
is the band structure which is derived from Eq.(2) and 
used to compare with our method. In Fig. (2a), we can 



FIG. 2: The numerical results of Fig.l. (a) is the band 
structure derived from Eq.(2) and in which the bold line 
is the W = 0.22^p line. (b) and (c) are the equal fre- 
quency contour line of propagation modes in k space and the 
min(| Imffc) 7^ 0) vs. k y figure, alternatively. The circle in 
(b) denotes the incident light of which u — 0.2^p. Both of 
them are derived from Eq.(6) when lu — 0.2^p, k z = and 
k v is scanned from — — to — . 

Besides, we can find that there are two propagation 
modes toward right when k y is a fixed number in Fig. (2b). 
These modes are similar to TE and TM modes in 2D 
isotropic PC, however, they can not be distinguished in 
3D PC, we just plot them directly. Furthermore, C\ sym- 
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metry exists in Fig. (2b) but not in the figure of real part 
of complex k x . The reason is the real number solutions 
of k x are the k x of the propagation modes which are the 
solutions of bulk system in which the C\ symmetry exist. 
However, the above is not correct when k x are complex 
numbers, because the complex means that there is an in- 
terface destroying the C\ symmetry and facing x direc- 
tion in the system as well. Therefore, all the evanescent 
modes of which k x are complex numbers just exist near 
the interface and their penetration depths correspond to 
2tt/ |Im(fc x ) ^ 0| owing to e* kr = e tk ' r e l ^ x , where R 
and / denote real and imaginary parts, alternatively. The 
most remarkable one of the complex k x relates to the 
longest penetration depth denoted as Xlpd {k y , k z ,co), 
because almost nothing but the propagation modes can 
exist in this system when the distance from the detecting 
position to the interface is larger than Xlpd {k yi k Zl uo). 
Therefore, a semi-infinite system can be treated as two in- 
dividual regions: surface and bulk regions, all the evanes- 
cent modes just exist in the surface region of which the 
width is As definded as max (Xlpd (k y , k z , u>o)), where 
loq is a fixed frequency. For a finite size PC, if the ef- 
fect of corner is not important, X$ decides the smallest 
size of PC. If the size is smaller than the smallest one, 
the system no longer can be treated as a periodic struc- 
tured media. Figure(2c) is the figure o{ci/Xlpd vs - k y at 
uj = 0.2^p. This figure indicates that the ci/Xlpd drops 
to zero quickly when k y is located at the edge of contour 
in Fig. (2b). This kind of situation arises while the state 
located at the edge of contour changes from propagation 
mode to evanescent mode. Besides, because \k y \ < 0.2^ 
when the incident light is a propagation mode in vac- 
uum, we can find that ci/Xlpd > 0.7. Therefore, the 
longest penetration depth is a/0.7 for all incident light 
perpendicular to z direction. 

In conclusion, because Eq.(6) is a k x eigenvalue equa- 
tion when ui, k y and k z are provided, the u> can be a real 
number at any time, and e and fi can be the function of 
ui, k y and k z or complex tensors. In addition, since most 
of k x are complex numbers, the minimum of |Im(fca;) ^ 
must exist, and this value will decide how large a PC is 
able to treated as a single crystal if the influence of corner 
is not important. Therefore, one of the issues we proceed 
to research is the influence of corner. We thank Prof. 
Ping Shcn for his opinion to excite us to find out the 3D 
formula EPWE method. 



I. APPENDIX 

The eij shown as below is the abbreviation of CG-G',ij- 




5i.il =G x e zz 1 (k + G% + (k + G) z e- y 1 (k + G') z - 

G x e zy \k + G') z -{k + G) z e xz 1 {k + G') y + uj 2 » yx , 
B h2 i =G x e- y \k + G') z + (k + G) y e x l(k + G') y - 

G x e yz \k + G') y -(k + G) v e xy L (k + G% + cu 2 [i zx , 
B 2 ,n =G x e' zx 1 {k + G') z + (k + G) z e xz G x - 

G x e- z l G' x - (k + G) z e~ l x (k + G% + c 2 Mw , 
B 2 ,i2 =G x e~lG' x + (k + G) z e~l(k + G') y - 

G x e~ zx \k + G') y - (k + G) z e~yG' x + lo 2 ^, 
B 2 ,2i ^G x e-lG' x + (k + G) y e x l{k + G%- 

G x e yx (k + G') z — (k + G) y e xz 1 G' x + Lo 2 ^ zy , 
-82,22 =G x e~ x (k + G') y + (k + G) y e X yG' x — 

G x e~lG' x -(k + G) y e-l{k + G') y + lu 2 ^ zz , 



C hll =e~J(k + G% - e- y \k + G') z , 
Ci,2i =e m 1 (k + G') z -e- z 1 (k + G') y , 



C2,ll 


=t Z x(k + 


G'M 


-(k + G) z e x l- 


e zz(Gx + G' x ), 


C2,12 


=e zy {Gx 


+ G' X ) 


-e- x \k + G') y 


-(k + G) z e x ^, 


C2,21 


=e yz (G x 


+ G' X ) 


-e yx \k + G') z 


-{k + G) y e xz 1 , 


C2,22 


= e yx (k + 


G') y ■+ 


- (k + G) y e xy - 


e yy(G x +G' X ), 



D 11 =(k + G) z e- z 1 -(k + G) y e zz \ 
D12 =(k + G) y e Z y -(k + G) z e~ y , 



E 1 =(k + G) y e zy 1 (k + G% + (k + G) z e- Z 1 (k + G') y - 

(k + G) y e- Z 1 (k + G') y -(k + G) z e yy (k + G% + lo 2 ^.. 

E 2 ,u =(k + G)ye zz 1 G' x + (k + G) z e yx \k + G') z - 

(k + G) y e- zx 1 {k + G') z - (k + G) z e~}G' x + uj 2 fi xy , 

E 2 ,i2 =(k + G) z e yy G' x + (k + G) y t- Z lik + G') y — 

(k + G) z e-l(k + G% -(k + G) y e~yG' x + lu 2 ^ xz . 
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